# Monte Carlo Integration
# Author: Ajay Rawat
#         Indira Gandhi Centre for Atomic Research 
# Date  : Thu Jun 9, 2011
import math
import numpy as np
import matplotlib.pyplot as plt
a = 0.
b = math.pi
Trials = 10001
f = lambda x: 0.5*x*math.sin(x)
value = []
for j in range (1, 101):
    print j
    xsum = 0.
    for i in range (1, Trials):
        x = (b-a)*np.random.rand() + a
        y = f(x)
        xsum += y
    ans =  (b-a)*xsum/(Trials - 1)
    value.append(ans)
    print 'Monte carlo Integration value:' 
    print value
plt.hist(value, bins = 10)
plt.show()
